function [s_xx,s_yy,s_xy] = ex33_AnalyticalStress(x,y,a)
% function [s_xx,s_yy,s_xy] = ex33_AnalyticalStress(x,y,a)
% This function computes the stress tensor = [ s_xx s_xy; s_xy s_yy]
% x,y = nodal coordinates
% a   = radius of the hole

eps   = 1.0e-7;
a     = abs(a);
r     = sqrt(x*x+y*y);

if abs(r-a)<eps
    r = a;
end

if abs(y)<eps
    theta = 0;
elseif abs(x)<eps
    theta = 90*pi/180;
else    
    theta = atan(y/x);
end

aux = a*a/r/r;

if r<a
    disp('r must be larger than a')
    s_xx = 0;
    s_yy = 0;
    s_xy = 0;
else
    s_xx = 1 - aux*(3/2*cos(2*theta) + cos(4*theta)) + 3/2*aux*aux*cos(4*theta);
    s_yy =   - aux*(1/2*cos(2*theta) - cos(4*theta)) - 3/2*aux*aux*cos(4*theta);
    s_xy =   - aux*(1/2*sin(2*theta) + sin(4*theta)) + 3/2*aux*aux*sin(4*theta);
end

